;+
; NAME:
;       IMIR_MOSAIC
;
; PURPOSE:
;       Generates a mosiac of the I_MIR fits for the BGPS sources
;       included in an rb3 input file.  This mosaic may be used to
;       derive a model I_MIR(l,b) based on the emission density
;       derived from the IRDC model fitting.
;
; CATEGORY:
;       IRDC Distance Modeling
;
; CALLING SEQUENCE:
;       IMIR_MOSAIC
;
; INPUTS:
;       NONE
;
; OPTIONAL INPUTS:
;       NONE
;
; KEYWORD PARAMETERS:
;       PLOT    -- Generate a PostScript plot of the mosaic.
;       VERBOSE -- Be verbose in output.
;
; OUTPUTS:
;       NONE
;
; OPTIONAL OUTPUTS:
;       SUCCESS -- Was the routine successful in doing its dirty deed?
;
; COMMON BLOCKS:
;
;
; MODIFICATION HISTORY:
;
;       Created:  04/09/12, TPEB -- Initial version.
;       Modified: 06/25/12, TPEB -- Codified version as a PRO able to
;                                   be called from another routine
;                                   (and changed name).
;
;-

;; Function to generate a model IMIR based on a radiation field fit
FUNCTION IMIR_MOSAIC_MODEL, im, P, mt, astr, hd
  
  COMMON FFORE_BLOCK,n,rb3,R0,d,R,Z,tau,f_data,sig_f,corr,farlist,do_Tfit,$
     lpstr,rho_hi,rho_h2,rho_star
  
  ;; Parse input
  P = double(P)
  A  = P[0]
  B  = P[1]
  C  = P[2]
  NG = P[3]
  TD = P[4]
  mt = P[5]
  GW = P[6]
  
  sz = size(im,/DIM)
  x = cmReplicate(lindgen(sz[0])+1, sz[1])
  y = TRANSPOSE(cmReplicate(lindgen(sz[1])+1, sz[0]))
  
  ;; Get galactic params
  defsysv, '!MW', exists = exists
  IF NOT exists THEN galactic_params 
  ;; Using !MW parameters -- Goes into COMMON block!
  R0 = !MW.R0
  d = dindgen(!MW.NBINS*5.)*!MW.BINSIZE + !MW.BINSTART
  
  model = im*0.
  
  ;; Find which pixels are nonzero!
  
  inz = WHERE(im NE 0., nnz)
  IF nnz EQ 0 THEN BEGIN
     message,'Input image is blank!',/cont
     RETURN,model
  ENDIF
  
  t_start = systime(1)
  ;; Need to loop over relevant pixels here...
  FOR k=0L, nnz-1 DO BEGIN
     
     IF (k+1) MOD 1000 EQ 0 THEN BEGIN
        total = (systime(1) - t_start) / (double(k+1) / double(nnz) )
        etd = systime(0, t_start + total)
        
        message,'Pixel #'+string(k+1,format="(I8)")+$
                ' of '+string(nnz,format="(I8)")+' ['+etd+']',/inf
        
     ENDIF
     
     xy2ad,x[inz[k]],y[inz[k]],astr,gl,gb
     omni_lbd2rz, gl, gb, d, r, z
     
     IF mt EQ 0 || mt EQ 2 THEN BEGIN
        rho_hi = hi_density(R,Z)
        rho_h2 = h2_density(R,Z)
        rho_star = star_density(R,Z)
     ENDIF
     
     ;; Select where to get model
     CASE mt OF
        0: junk = BGPS_FFORE_MODEL(1.d3, A, B, C, NG, rho)
        1: junk = BGPS_FFORE_SIMPLE(1.d3, A, B, C, rho)
        2: junk = BGPS_FFORE_MODEL2(1.d3, A, B, C, rho)
        3: junk = BGPS_FFORE_SIMPLE2(1.d3, A, B, C, NG, GW, rho)
        ELSE:  message,'Model Type '+string(mt,format="(I0)")+' not valid.'
     ENDCASE
     
     inten = total(rho)
     model[inz[k]] = inten
  ENDFOR
  
  ;; Normalize to get close...
  model = model / max(model,/NAN) * max(im,/NAN)
  
  undefine,x,y
  RETURN,model
END


;;==================================================================
;; MAIN ROUTINE
;;==================================================================
PRO IMIR_MOSAIC, SUBSET=subset, PLOT=plot, VERBOSE=verbose, SUCCESS=success
  
  ;; Parse keywords
  IF n_elements(subset) EQ 0 THEN subset = ''
  IF KEYWORD_SET( verbose ) THEN verbose=1b ELSE verbose=0b
  success = 0b
  
  ;; Restore the desired RB3 file for use here...
  rfn = './irdc_dist_model/bgps_rb3'+subset+'.sav'
  IF ~ FILE_TEST(rfn) THEN BEGIN
     message,'File '+rfn+' does not exist.  Returning...',/cont
     RETURN
  ENDIF
  restore,rfn,VERBOSE=verbose
  n = n_elements(rb3)
  
  
  ;; Create the blank image for the Galactic Plane, header & astrometry
  im = fltarr( fix(45./(7.2/3600.))+2, fix(2.5/(7.2/3600.))+2)
  
  mkhdr, hd, im
  bgsz = long(size(im,/DIM))
  make_astr,ast,CD=[[-0.00199999986216d,0.],[0.,0.00199999986216d]],$
            CTYPE=['GLON-CAR','GLAT-CAR'],CRPIX=bgsz/2L,CRVAL=[42.5d,0.d]
  putast,hd,ast
  
  ;; Loop through sources in RB3 and add their I_MIR fit postage stamp
  ;; to the mosaic
  t_start = systime(1)
  FOR i=0L, n-1 DO BEGIN
     
     ;; Read in postage stamp and get relevant astrometry / pixel
     ;;    information
     fn = './local/ubc_glimpse_proc/Imir'+string(rb3[i].cnum,format="(I04)")+$
          '.fits'
     imir = readfits(fn,ihd,/silent)
     sz   = size(imir,/DIM)
     xy2xy,ihd,hd,0,0,ix,iy
     ix = long(round(ix))
     iy = long(round(iy))
     
     ;; Add postage stamp!
     im[ix>0 : ((ix+sz[0])<bgsz[0])-1, iy>0 : ((iy+sz[1])<bgsz[1])-1] = imir
     
  ENDFOR
  IF verbose THEN $
     print,'Mosaic creation time (seconds): ',systime(1) - t_start
  
  ;; Write out the mosaic to a FITS file
  ;; Add Comments AND/OR History?
  
  IF verbose THEN print,hd
  filename = 'Imir_mosaic'+subset+'.fits'
  sxaddpar,hd,'FILENAME',filename,BEFORE='CTYPE1'
  sxaddhist,'I_MIR Mosaic image produced by IMIR_MOSAIC.pro for BGPS '+$
            'IRDC distance modeling',hd
  
  writefits,'./local/output/'+filename,im,hd
  success = 1b
  
  
  ;; Make the crazy plots, if desired.
  IF KEYWORD_SET( plot ) THEN BEGIN
     ct = 13
     colors = [42,84,127,170,212,255]
     rev=0b
     
     myps,'./irdc_dist_model/analysis_plots/glon_vs_Imir.eps',ct=ct,xsize=15,$
          REVERSE=rev,ysize=10
     
     multiplot,[1,4],gap=0.02
     cs = 1.0
     
     ind = WHERE(im NE 0., nind)
     range = alog([10,max(im[ind])])
     plotmap,im,hd,range=range,/log,charsize=cs,ct=ct,REVERSE=rev,xr=[65,42.5]
     
     print,exp(range)
     cgLoadct,ct,/silent,REVERSE=rev
     cgColorbar,range=exp(range),title='I!dMIR!n GLIMPSE Intensity [MJy/sr]',$
                charsize=cs,/xlog,xtickformat="(F0.1)",divisions=6,$
                position=[0.15,0.95,0.90,0.97]
     
     multiplot,/doxaxis,/doyaxis
     
     plotsym,0,0.5,/fill
     plot,rb3.l,rb3.i_mir,xr=[65,42.5],psym=8,xtit='Galactic Longitude [deg]',$
          ytit='I!dMIR!n [MJy/sr]',/xst,/nodata,charsize=cs,$
          color=cgColor('opposite'),yr=[0,200]
     
     lev = exp((dindgen(6)+1.d)/6.*(range[1]-range[0]) + range[0])
     print,lev,format="(F0.1)"
     
     FOR jj=0,5 DO vline,/horiz,lev[jj],linestyle=0,thick=3,color=colors[jj]
     
     oplot,rb3.l,rb3.i_mir,psym=8,color=cgColor('opposite')
     
     multiplot,/doxaxis,/doyaxis
     
     plotmap,im,hd,range=range,/log,charsize=cs,ct=ct,REVERSE=rev,xr=[42.5,20]
     
     multiplot,/doxaxis,/doyaxis
     
     plotsym,0,0.5,/fill
     plot,rb3.l,rb3.i_mir,xr=[42.5,20],psym=8,xtit='Galactic Longitude [deg]',$
          ytit='I!dMIR!n [MJy/sr]',/xst,/nodata,charsize=cs,$
          color=cgColor('opposite'),yr=[0,200]
     
     lev = exp((dindgen(6)+1.d)/6.*(range[1]-range[0]) + range[0])
     print,lev,format="(F0.1)"
     
     FOR jj=0,5 DO vline,/horiz,lev[jj],linestyle=0,thick=3,color=colors[jj]
     
     oplot,rb3.l,rb3.i_mir,psym=8,color=cgColor('opposite')
     
     myps,/done,/mp
     
  ENDIF
  
END
